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Abstract 

A radial basis function (RBF) method based on matrix-valued kernels is presented and an¬ 
alyzed for computing two types of vector decompositions on bounded domains: one where the 
normal component of the divergence-free part of the field is specified on the boundary, and one 
where the tangential component of the curl-free part of the field specified. These two decom¬ 
positions can then be combined to obtain a full Helmholtz-Hodge decomposition of the field, 
i.e. the sum of divergence-free, curl-free, and harmonic fields. All decompositions are computed 
from samples of the field at (possibly scattered) nodes over the domain, and all boundary con¬ 
ditions are imposed on the vector fields, not their potentials, distinguishing this technique from 
many current methods. Sobolev-type error estimates for the various decompositions are pro¬ 
vided and demonstrated with numerical examples. Radial Basis Functions; Kernel Methods; 
Vector Decomposition; Divergence-free Approximation; Curl-free Approximation. 


1 Introduction 

In the literature the phrases “Helmholtz decomposition,” “Hodge decomposition,” and “Helmoltz- 
Hodge decomposition” are used to describe a variety of vector decompositions in which a given field 
f is written as a sum of divergence-free and curl-free fields. We will refer to any such decomposi¬ 
tion as a Helmholz-Hodge decomposition (HHD). These decompositions are fundamental to many 
applications, from fluid dynamics and electromagnetics, to computer graphics and imaging. Each 
component plays an essential role in the underlying application. For example, the incompressible 
Navier-Stokes’ equations describe the dynamics of an incompressible fluid, the velocity field of the 
fluid is divergence-free while the (hydrostatic) pressure is curl-free. This fact is exploited in projection 
methods, which are the dominant strategy employed for numerically solving these equations [511291. 
A more general version of such a decomposition is given by the Hodge Theorem m , which implies 
that vector fields f on a compact domain fi C can be split into the sum f = w + Vp + Vh, where 
w is divergence-free and tangent to the boundary, Vp is curl-free and normal to the boundary, and 
the scalar function h is harmonic. This “full” HHD is used in graphics for detecting singularities 
(e.g. sinks, sources, and vortices) in vector fields that arise in various disciplines [25] , 

Several techniques exist to compute HHDs, with most making use of the vector field sampled 
on a mesh or grid. The standard approach employed is to recast the problem in terms of a Poisson 
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equation for a potential function p. More specifically, given a vector field f, one numerically solves 
A p = V • f, using, for example, finite difference or finite element methods. It follows then that f 
is the sum of Vp (which is curl-free) and f — Vp (which is approximately divergence free). One 
drawback of this approach is that in many applications it is not clear how to impose the correct 
boundary conditions on the Poisson problem for the potential p. This is in part because the bound¬ 
ary conditions are typically imposed on the divergence-free or curl-free fields directly, not on the 
potentials for these fields. For example, with regard to solving the incompressible Navier-Stokes 
equation, standard projection methods require a decomposition by calculating a pressure p as the 
solution of a Poisson problem. However, the pressure does not have a boundary condition as it plays 
the role of a Lagrange multiplier, with its value being whatever it has to be to make the velocity 
field divergence-free [5] - 

Other techniques for decomposing vector fields use basis functions that are customized to split 
into analytically divergence- and curl-free parts. These methods avoid having to explicitly solve a 
Poisson problem, but do require solving some other type of problem (e.g. an interpolation problem). 
Examples on periodic domains include those utilizing wavelets [9] , and meshless kernel methods such 
as spherical basis functions mm For domains with boundaries, a meshless radial basis function 
(RBF) method was developed for numerically solving certain static fluid problems (see [26 , 31]), 
with a by-product of this approach being a method for computing a certain type of decomposition. 

In this paper we develop and provide error estimates for a meshless RBF method for computing 
two standard vector decompositions on bounded domains in R d - 2 : one where the normal component 
of the divergence-free part of the field specified on the boundary, and one where the tangential 
component of the curl-free part of the field is specified. These decompositions can then be combined 
to compute the full HHD on a bounded domain. Our approach utilizes matrix-valued RBFs that 
split into analytically divergence-free and curl-free parts. Each decomposition is obtained by solving 
a generalized interpolation problem, with the boundary conditions appearing on the velocity field 
variables and not on the potentials, and gives rise to a positive definite linear system of equations. 
While we never work with the (vector and scalar) potentials of the components of the decomposed 
field directly, these potentials can be easily recovered at no added computational cost. Our method 
provides accurate decompositions, but does require global information. As such, a drawback, as 
is the case with many global kernel-based methods, is expense. We hope this can be mitigated 
by employing approaches similar to those in the scalar kernel theory, such as using a multiscale 
approach [TO] or by employing a localized basis El HE but this will be reported on separately. 

As noted above the technique described in [26] SI] also gives rise to methods for computing certain 
vector decompositions in R d . In fact, a vector decomposition as in Proposition 1 was obtained in 
EE In these papers the authors use “combined kernels”, which are constructed by incorporating 
a d x d divergence-free kernel with a scalar RBF to obtain a larger (d + 1) x (d + 1) kernel. Our 
approach is different in that instead of combining kernels to make a larger one, we sum kernels with 
properties to match the HHD, which results in a diagonal d x d matrix-valued kernel. Though not 
obvious at first appearance, it can be shown that the techniques are in fact equivalent for a certain 
choice of the scalar kernel in the combined method. However, we approach the problem from a 
different perspective—instead of using a combined kernel that sets out to model the components of 
the vector field with separate kernels, we model the field directly with a single kernel that splits 
naturally. A practical by-product of this approach is that a large portion of the interpolation matrix 
becomes block-diagonal, which gives savings in terms of storage and computational efficiency. Where 
there is overlap in our work with previous work, we offer improvements in error estimates in terms 
of the order of approximatioi)]] and the domains on which they apply. We also include a vector 
decomposition not treated with kernel methods before (as described in Proposition 2) and develop 
the first kernel method for computing the full HHD. 

The paper is organized as follows. Section [2] contains the necessary preliminaries on function 

1 Previous work derived estimates measured in the H 1 norm. We extend this to L 2 , which gives an extra order of 

approximation. 
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spaces and vector decompositions. In Section[3]we give background information on scalar and matrix¬ 
valued RBFs. Next, the construction of our kernel decompositions are described in detail in Section 
[2 Error estimates and numerical experiments are presented in Sections [5] and [G] respectively. We end 
the paper with some concluding remarks regarding decompositions with other boundary conditions. 


2 Preliminaries 

We will distinguish between scalar and vector valued functions by denoting the latter in bold-face. 
We denote the gradient and divergence in the usual way, i.e. V and V-. The curl operator on three 
dimensional fields will be denoted by curl(f). Given a scalar valued function / : R 2 — > R, we will 
use the same notation for curl(/) := (—d y f, d x f ) — this should cause no confusion. We will let H 
denote a connected open domain in R d with boundary T of Holder class C" 1,1 for some nonnegative 
integer m. 

2.1 Function spaces 

The function spaces we will work with are all Hilbert spaces: L 2 {^L) will denote the space of square 
integrable functions on H, and will denote the space of all vector fields with in L 2 (VL). Given 

s > 0, we let H s (£l) denote the Sobolev class of functions on H with smoothness s, and denote its 
vectorial analogue by H s (fi). When the underlying domain is R d , we use the Fourier transform form 
of the inner product in these spaces. For example, the inner product on H s (R d ) is given by 

(F §) H s (M d ) := f IM(i + M 2 ) s dw, (l) 

JR* 

where f denotes the Fourier transform of f and |w| denotes the Euclidean length of w € R d . We will 
also need the space of functions H* s (R d ), which is endowed with the inner product 

f ~— fl + |w | 2 ) s+1 

(/.ffW)-/ 2 du>. (2) 

Js. d M 

It can be shown that H s (M. d ) is a subspace of H s (R d ) and that ||/||iy»(R<*) < ||/||^»( R d) for all 

/ £ H s (M. d ) [T5], Proposition 2]. The space H s (R d ) is defined in an analogous way. 

We denote the L 2 (r) inner product by (-, •). Sobolev spaces on the boundary T can be defined in 
various ways. If the boundary is C 171 ’ 1 , then to define H S {T) with 0 < s < m + 1 one can use charts 
and a partition of unity (see, for example [191 Section 1.3.3]). For s > 0, we let H~ S {T) denote the 
dual space to H S {T), and the vector-valued cases for these spaces will be denoted in bold-face. 

Our arguments later will require standard operator interpolation on Sobolev spaces. A concise 
treatment of what we need can be found in [U Ch. 14]. For the interpolation arguments on boundary 
spaces, we will use the following fact from pii . Theorem 7.7]: Let 0 < 6 < 1. For all si, S 2 £ R with 
si > S 2 we have 

[H Sl (T),H S2 (T)] e = #( i-0)*i+e*2( r ) j ( 3 ) 

with equivalent norms, where [H Sl (T), H S2 (r)]g is the interpolation space with parameter 8 between 
H S1 (T) and H S2 {T). 

Lastly, we will make use of the following norms, which are both equivalent to || • ||h s ( 0 ) f° r a ll 
s > 1 when T is at least 

11 u 11 l 2 (£ 2 ) + || cu rl( u )||H'>- 1 (fi) + ||V • u||^-a-i(Q) + || u ■ n||^ s _ 1/2(r ), 

II u IIl 2 ( 0 ) + ||cur 1 (u)||^ s _ 1(Q ) + IIV • u||^ s _ 1(n) + ||u x n||^_ 1/2(r) . 
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(4) 

(5) 





For integer s, see p~Bl Corollary 3.7, pg 56] for © and [71 Proposition 6’, pg. 237] and the proceeding 
remarks for ©• The fractional cases follow from standard interpolation arguments. Though stated 
here for d = 3, similar results hold in the two dimensional case. 


2.2 Vector Decompositions 

The Helmholtz-Hodge decomposition for vector fields in L 2 (R d ) can be easily described in terms 
of the Fourier transform. A field f £ L 2 (M d ) is divergence-free if and only if up^tiuf) = 0 almost 
everywhere, and f is curl-free if and only if f(w) = ioh(uj) for some h £ H 1 ( R d ). Letting T~ x : 
L 2 (R d ) —> L 2 (R d ) denote the inverse Fourier transform, the operators 


D f_77—1 

r div L • J 



Pcurl? '■= F 1 



( 6 ) 


are projections on L 2 (R d ), with Pdiv f divergence-free, P cur /f curl-free, and P^f T P cu wf- With 
this, f = Pdiv f + Pcurtf uniquely decomposes f into L 2 (R d )-orthogonal divergence-free and curl-free 
fields. Further, P*„ and P cur i are also orthogonal projections on any space whose inner product is 
of the form 

(f,g)*= [ g (w)ip{ui)du}, (7) 

J K d 

where the weight function <p > 0 is measurable—this includes all Sobolev spaces H s (R d ) and H s (R d ). 
For fields on bounded domains we will focus on the two fundamental decompositions given in the 
following propositions. 


Proposition 1. Let tl C R d be a connected Lipschitz domain, f £ L 2 (f2) be such that V f £ L 2 (f2), 
and let g £ H - 1 ' 2 (r) satisfy ( g , 1) = 0. Then one has the unique decomposition f = w + Vp, where 
p £ H 1 (S7), and w £ L 2 (f2) satisfies V • w = 0 with w • n = g on T. The function p is uniquely 
determined up to a constant, and satisfies the bound 


biffin) = ||Vp|| i2 (n) < C (||V • f||i 2 (n) + ||f • n - g\\ H -i/ 2 ^ r) ) , (8) 

where C is some constant independent oft. When g = 0, w and Vp are orthogonal in L 2 (fi). 

Proof. Since the divergence of f is in L 2 (f2), f has a well-defined normal boundary component 
f ■ n £ H~ l / 2 (T) satisfying Green’s formula (see [TB] Theorem 2.5]). Thus we can consider the 
following weak Neumann problem 

(Vp, Vr) = (—V • f, v) + (f • n — g, v) V v £ H 1 (tl). 


Standard Lax-Milgram theory dictates that the solution p is continuous with respect to the data, 
giving © ( see, for example, PI Proposition 1.2]). The field w := f — Vp has the other properties 
listed above. □ 


An important by-product of this decomposition in the case g = 0 is the Leray projector Pl and its 
orthogonal complement P^, defined by Pl f := w and Pft := X7p. 

The next decomposition splits a vector field into a divergence-free field and a gradient field 
normal to the boundary. Note that Vp is normal to the boundary if and only if p|r is constant on 
each of the connected components of T, which we denote by rn,ri,...,T/<-. The following is from 
Corollary 5' in © pg 224], 

Proposition 2. Every f £ L 2 (fl) admits the unique orthogonal decomposition f = w + Vp, where 
p £ = {u £ H 1 (LI), = constant , * = 0,..., K}. The vector field w is divergence-free and 

perpendicular to S7p in L 2 (f2). 
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2.2.1 Potential Functions and Extensions 

In what follows we require w (the divergence-free term of f) to be expressed as w = curl(i/>) in 
the case of d = 3 dimensions (or w = curl(t/i) when d = 2)0 We will also need a well-defined 
continuous assigment w —> ip. This requires some mild assumptions on f l in the event that f 1 is 
multiply connected. Specifically, we assume that fl can be made simply connected by a series of 
non-intersecting “cuts” £i,..., £„, where Ej C is a smooth variety (see for example m pg. 2i7]). 
On such an fl, we have the following: 

Propositions. A given w £ L 2 (fl) is an element o/curl(H 1 (fl)) if and only if w satisfies Vw = 0 
and f r w ■ ndT = 0 for all i = 0 ...K. Of all possible potential functions, there is a unique 
ip £ H 1 (fl) such that w = curl (ip) satisfying 

V • if = 0, if ■ n = 0, (ip ■ n, l)s 4 = 0, i = 1,..., n. (9) 

Finally, we have the bound ||'0||H 1 (n) < C1 |w||l 2 (q) for some C independent of w. 

Proof. The first claim is Corollary 4 from [7] pg. 224], and the unique assignment follows from 
Remark 4 proceeding the corollary. For continuity, note that curl(H 1 (Sl)) endowed with the L 2 (fl) 
norm is closed [7] pg. 222, Proposition 3]. Now let V denote the subspace of fields ip £ L 2 (fl) 
satisfying ©. By [7J pg. 225, Proposition 4], V is closed in L 2 (fl), so V flH 1 (fl) is closed in H 1 (fl). 
Using this one can show that the operator T : curl(H 1 (fl)) —> V fl H 1 (U) given by T w := ip is a 
closed map, and therefore continuous. □ 

This leads to potential functions for our decompositions that satisfy the following regularity 
result. 

Proposition 4. Let r be such that 0 < r < m and let f £ H T (fl). Then the decompositions 
in Propositions [7] and[M can be written as f = curl(i/>) + Vp, for uniquely determined potentials 
p £ ff T+1 (fl) and ip £ H r+1 (fl). For the decomposition in Proposition [7] with g £ i7 T_1 / 2 (T) 
satisfying ( g , l)r 4 = 0 on each connected component ofT, these potentials satisfy 

Ibll^+qn) < C'(||f|| H ’-(n) + Hflili?’—va(r))) ||V ; IIh t + 1 < C(l|f||H-(n) + lbllfl^-i/ 2 (r))j (1°) 

Similar bounds (with g = 0) hold for the decomposition in Proposition [H 

Proof. Let r be a nonnegative integer. In the case of Proposition []] with g = 0, existence and 
uniqueness of ip follows from [7] page 224, Corollary 5] and the proceeding remarks. The Proposition 
[2] case follows from [7] page 224, Corollary 5']. The additional regularity of the boundary gives 
regularity of these potentials (see, for example [TJ page 236, Corollary 7]). Recall that V denotes 
the subspace of fields ip £ L 2 (U) satisfying ©, and V is closed in L 2 (U), so V fl H T+1 (fl) is 
closed in H T+1 (f2). From this one can show that the assignment f —> ip is a well-defined closed 
map, and thus obtain the bound for ip in m- The scalar potential p is unique if we require 
J n pdx = 0. In a similar fashion as above, the bound for p follows from the fact that the space 
H T+ \n) fl {p £ L 2 (fl) I J Q pdx = 0} is closed in iJ r+1 (U). The fractional cases can be handled 
using standard interpolation arguments. 

To handle the case g ^ 0 from Proposition [TJ let p g be the solution of the problem 

-A p g = 0 in U, ^ = -g on T, 

Note that that w g := — S7p g is divergence free. Since w g is divergence-free and w g ■ n = g satisfies 
the conditions in Proposition [3l w ff = curl(^ g ) for a unique ip . Letting f = curl(t/> 0 ) + V(po) 

2 Since our results will hold in two and three dimensions, throughout the remainder of the paper we will concentrate 
specifically on more complicated the d = 3 case to avoid constantly distinguishing between these two cases. 
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denote the decomposition of f from Proposition [T| with g = 0, where the potentials are the unique 
potentials from above satisfying (fTOl) with g = 0, the desired potentials are given by ip := if ) 0 + ip g 
and p := p 0 +p g . 

The bound m will follow from bounding ip g and p g . Since g £ H T 1 // 2 (r) and the domain is 
assumed smooth enough, we get the regularity bound [18l Theorem 1.10] 

||w ff || H -(fi) = ||Vp g || H -(fi) < lb 9 ||i/-+i(n) < C\\g\\ HT -i/2(j~y 

Using this with Proposition [3] ip g satisfies the bound ||V’ 9 ||H 1 (fi) < C||w 9 || L2 (n) < C\\g\\ H -i/ 2 (Yy 
For higher regularity, we use (U) with s = r + 1 to finish the proof: 

HV’glllb+^n) ~ lll^sllln — G (ll W sllH 1 (n) + ll w sllH T (n)) ^ ( ^ , |l w sllH T (f2) ^ C , ||fl , llHr-i/ 2 (r r 

□ 

We remark that the existence of these potentials is only used for theoretical purposes. The choice of 
cuts and the conditions (O plays no role in implementing the kernel-based decomposition presented 
later. However, potential functions for each term in the kernel decomposition will be readily available. 
Next we use these potentials to define an extension operator, which will be useful later. 

Lemma 1. Let g £ H T ~ 1 ^ 2 (T) satisfy (p, l)r 4 = 0 on each connected component of T, and let f = 
w + Vp denote the corresponding vector decomposition from Proposition QJ Given LI C R d satisfying 
the assumptions preceeding Proposition [7J there exists an extension operator E : H r (U) —> H r (R d ), 
for all t satisfying 0 < r < m, such that 

Ei |n = f, PdivEi\n = w and P cur iEi\n = Vp, (11) 

and is continuous in the sense that ||i?f ||g x ^ Rd ^ < C (||f||H x ( 0 ) + lblli/ T - 1 / 2 (r)) • 

Proof. Let p and ip denote the unique potentials for a given f £ H T (U) in Proposition!!] These can be 
extended using Stein’s continuous extension <£ : H T+1 (Ll) —> H r+1 { R d ), which we note is universal in 
the sense that <£ does not depend on r [25} Chapter 4]. We will interpret € : H r+1 (f2) —> H r+1 (R d ) 
as € applied component-wise. We can then define the extension Ef := curl(£t/>) + V(£p, which 
satisfies m■ Lastly, (usd gives us that E is continuous: 

ll^ f llg T(R d) = (b x + |w£p| 2 ) ( 1 + |^j 2 ' > - duj 

< j (|^| 2 + I^P| 2 ) (1 + | W | 2 ) T + 1 du) = ||(£^||Hr + l( R d) + ||£p||ffT + l(H<<) 

< CII’AllH^+qn) + ^ ^ (lbllH T (n) + llffll^-i/2(r)) ■ 

□ 

These same arguments can be repeated to establish a continuous extension satisfying m for the 
decomposition in Proposition [2] 

3 Radial Basis Functions and Related Kernels 

A kernel <p : R d x R d —> R is positive definite if given any finite set of unique points X = 
{xi, X 2 ,.. ■, Xn} C R d , the associated Gram matrix with entries Aij = cp(xi,Xj) is positive defi¬ 
nite. The typical Ansatz for interpolation of function / over the points X with such a kernel is to 
find an interpolant of the form 

N 

S f = *52 <!>(•, X j) C 3, ( 12 ) 

J=1 
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where the coefficients Cj are chosen so that s/1 Y = f\ x - Positive definiteness of the kernel ensures 
existence and uniqueness of the interpolant. If <J) is radial in the sense that 4>(x,y) = y>{\x — y |) for 
some univariate <p, then 0 is a radial basis function (RBF). It is common to simply write (j>{x,y) = 
(j)(\x — y |). Good references on RBFs are, for example, [5] fill [50]. 

For vector-valued approximations, there are matrix-valued kernels $ : R d x R d —> R d x R d . 
Interpolants to a vector field f : R d — > R d sampled at distinct points X = {x\,X 2 , ■ ■ ■ ,Xn} C R d can 
be constructed from these kernels as follows: 


N 


s f = 

3 =1 


3 > 


(13) 


where the vector coefficients c j £ R d are chosen so that Sf | Y = f| y . This leads to the following 
Nd x Nd linear system of equations: 


(14) 


$(xi,Xi) 

$(xi,xjv) 


Cl 


■fr 

$(aqv,xi) 

••• <F(XjV,Xjv)_ 


CJV. 


fjV 


A 


c 


f 


We say that is positive definite if the Gram matrix A in (fTdl) is positive definite for any distinct 
set of points X. It will be useful later to express this property in a block-style quadratic form. Since 
A is positive definite, we have 


^ Ck$(xk,Xj)cj = c T Ac > 0, (15) 

j,k 

with equality occurring if and only if Cj = 0, j = 1,...,N. 

Customized matrix-valued kernels leading to divergence-free and curl-free approximations were 
introduced independently by several researchers in the 1990s: mmm- In all cases the construction 
of the customized kernel is fairly simple. For example, letting <j> be an RBF on R 3 , we define 

$d.iv{x, y) = curl* curl y (cj>( \x - y|)I) and $ C url(x, y) = V^Vy (<j)(\x - j/|)I) , (16) 

where I is the 3-by-3 identity matrix, the subscript in the differential operators indicate which 
argument they act on, and the curl of a matrix is interpretted as having the curl operator act on 
the matrix column-wise. Note that X y cf) = —\7 x <j), so this simplifies to a form that readily generalizes 
to any M d : 

$div(x,y) ■■= (-AI + \/X T )cj)(\x - y\) and $ CU ri(x, y) := -VV T ^(|x - y\), 

where the differential operators act on x. It is easy to check that the second argument acts as a 
shift, e.g. &div(x, y) = &div{x — y). If </> is positive definite, 4and & C uri are both positive definite 
(see, for example mm)- Further, the kernel given by 

^ div "b *&curl = XcJA (17) 

is also positive definite because it is the sum of positive definite kernels. $ decomposes naturally 
into its divergence-free and curl-free components. Indeed, given Xj , c, £ R d , the identitieU 

®div(w) = (M 2 I - ulo t ) (j){(jj) and § C uri{oj) = (u juj t ) cj)(uj) 
imply that Xj)Cj mid ^curi^'t 

3 Here (j) denotes the d -variate Fourier tranform of the single argument function </>(| • |). 
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3.1 The Native Space 

From here on out, we let $ denote the matrix-valued kernel from O- Each positive definite matrix¬ 
valued kernel gives rise to a canonical reproducing kernel Hilbert space, commonly referred to as 
the native space for that kernel. The native space for $ is denoted by AA precise definition 
for A$(R d ) is not warranted here and we refer the interested reader to m Section 3]. $ serves as 
a reproducing kernel in the sense that if f is a vector field in Aand b £ K d , then 

(f^(-A)b)Af 4 (i^)= brf (a;) \/x£R d , (18) 

where (•, • )_y # (gd) denotes the inner product on A/$(R d ). 

It can be shown that if <f> £ C 2 (R d ) with Acf £ Li(M. d ), then the inner product in A/$(K d ) is 


(fi g)v*(R d ) 



■T 


fM §M 

M 2 0(w) 


du, 


(19) 


where f is the Fourier tranform of f and A/$(K d ) C L 2 (K d ) is identified with all functions finite 
in the associated norm (see (T51 Section 3.1]). It immediately follows that if the RBF <fi satisfies 
</)(u>) < C( 1 + M|) _T_1 for some constant C, then A$(M d ) is continuously embedded in H r (M d ). If 
in addition 

£(w) ~ (1 + M2)—\ (20) 

then A/$(M d ) = H r (R d ) with equivalent norms. 


3.2 Generalized Interpolation 

The reproducing kernel Hilbert space structure of the native space makes it possible to interpolate 
using a wide variety of continuous linear functionals. A concise treatment of this is given for scalar¬ 
valued RBFs in m Chapter 16], and generalizes in a straightforward way to the matrix-valued case. 
We summarize the main results we need below. 

Let A C Af<s>(Mr)* be a finite linearly independent collection of linear functionals, where A/$ (M®)* 
denotes the dual space to A/$(R d ). Given the data (A(f) | A £ A}, where f £ A/$(M d ), we look for a 
generalized interpolant to f of the form 


Sf = ^2 v A a A , 

AeA 

where ot\ £ R and each v A is the Riesz representer for A. The interpolation conditions A(sf) = A(f) 
V A £ A lead to a linear system, and as long as the functionals are linearly independent the problem 
is uniquely solvable. Further, Sf is perpendicular to f — Sf in A$(K d ), which gives us the following: 

ll f - s f|lA4(R d ) < IlfllAAfR 11 )) ll s f||jV*(Rrf) < ||f||v*(R'i)- (21) 

Note that since is a reproducing kernel for A$(K £i ), the Riesz representer for A can be written in 
terms of d>. For example, (fT51) shows that the evaluation functional defined by A(/) = b T f(xj) is 
represented in the native space as d>(-,a:j)b. Next we consider functionals involving Pdiv 

Proposition 5. Let x,n £ and define the functional f) := n T Pdi V i(x). Then v is continuous 
on A$(R rf ) and has Riesz representer ( hd^(-, x)n. 

Proof. First note that by (Hhd and (0, Pdi V is a projection on A$(R d ). Using this and the repro¬ 
ducing kernel property of $ we have 

Kf)l = 


\{Pdivf, ^(',a:)nj)jv*(K‘ I )l < ll^(£2;) n llAq.(R' 1 )ll-Pd™f||Aq,(R d ) < C||f||Aq.(R' i )- 




This gives us continuity. To verify the form of the representer, first note that the Fourier transform 
of g := $ div (-,x)n is given by 


g(w) = (M 2 I - w r )^)e“ T “n. 


Using this and (fl9l) . we have 


(f>g)A/i(R d ) 




f(w)e“ “dw = n T 


Pdiv f(w)e“ TtJ du = n 1 P div f(x). 


□ 


4 Kernel-based Decompositions 

In this section we show how to construct a kernel-based approximation to the decompositions dis¬ 
cussed earlier. We will also show how one easily obtains potential functions from the kernel approx¬ 
imation. 

4.1 Kernel Approximation with Divergence-free Boundary Conditions 

Given a target f on U and boundary target g , it is our aim to construct a kernel approximation s£ 
such that PdivSf and P CU riSf, which we can compute analytically, approximate the appropriate terms 
of the decomposition in Proposition [10 We will construct our kernel-based vector decomposition by 
requiring full interpolation on nodes X = {x\,X 2 ,... ,Xn} C U, while at the same time enforcing 
boundary conditions at a dense set of nodes Y = { 2 / 1 , J/ 2 , ■ ■ ■ ,2/m} C T. Although no repetition is 
allowed within each node set, X and Y can have a nonempty intersection. 

Letting e* £ denote the vector whose only nonzero entry is a 1 in the i th position, the 
interpolation functionals are given by A^(f) := for 1 < i < d, Xj £ X. The boundary 

functionals are given by v 3 (f) := Pd iv f (y : j), y 3 £ Y, where n y £ is the outward normal vector 
at y £ T. This gives a total of dN + M conditions to be met. The basis functions to be used are 
the Riesz representers of these functionals, which from the previous section are given by $(-,:Ej)ej 
and ^ d i V (-,yj)n yj , respectively. 

Using these as basis functions, our RBF approximation will take the form 

N d M N M 

Sf = ^ ^ ~b yj)n yj dj = ^^ ^(-,Xj)Cj + } ' < $>div{'iyj) n yjdj, (22) 

i=i »=i j=1 j =1 3 =1 

where the coefficents Cij, 1 < i < d have been consolidated into the vector unknowns c j for each 
j, as in (fT3|) . Letting f\x denote the dN x 1 vector whose j th d x 1 block is given by f (xj ), the 
interpolation conditions 1 and 2 above lead a linear system of the form 


A 

B ' 


c 


' fix ' 

B t 

C 


d 


. 9 \y _ 


where A is the matrix given in (114L B is given by 


B = 


^divix 1,2/i)n yi 


$div(xN,yi)n yi 


*&div (^ 1 , VM 
^div {%N 1 yM )Uy M 


4 We use the superscript t because when g = 0 the divergence-free portion is tangential to F. 
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and C is an Mx M matrix given by Cij = n^.^diviVi, Vj) n yj • Note that due to the diagonal structure 
of the kernel $ = A <f>I, the matrix A can be rearranged to be block-diagonal, with d identical N x N 
blocks along the diagonal. This not only reduces the cost of storing the interpolation matrix, but 
also makes it possible to solve (1231) using a more efficient Schur complement method than if the 
matrix A was dense [3j. 

Note that the interpolation matrix in (1231) is symmetric, and since we have taken the symmet¬ 
ric approach for generalized interpolation, it is also positive definite (and hence invertible) if the 
functionals involved are linearly independent [301 Section 16.1]. 

Lemma 2. The functionals in A = {A^ | Xj € X, 1 < i < d}A{vj \ yj G Y} are linearly independent. 

Proof. Suppose that some linear combination of the functionals in A sums to zero. This is equivalent 
to its Riesz representer vanishing, i.e. 

N M 

g := $(•, Xj)cj + ^2^div{-:yi)di = 0, 

j=i i=i 

where d; = nidi for some scalars di. Since the terms in the decomposition g = Pdi V g + Pcurig are 
orthogonal in Awe have ||P cur ig||^^ Rd j = 0. We also have 

||-PcuWg||jV*(R' J ) = ^ 'f di> C url{' 1 x j) c jj *&curl(', fffc) c fc)w$(R rf ) • 
j,k 

Using the native space inner product (1191) with the Fourier identities 

$c url(-,Xj)cj = (ww T )c^(w)e“^, $(-,x k )c k = c fc |w| 2 0(u;)e“* w , 

it follows that 

curl (’ 5 Ej )? ^curl (' ? %k (R d ) curl (* ? %j )c i ,$(-,x fc )c fc V 4>(R d ) . 

Thus the reproducing property of $ gives us 

||-fcwig||jV # (]R<i) = %j) c j %k) c k) AT$(R d ) = ^ " c fc *&curl ( x k i x j) c j ) 

j,k j,k 

and since & CU ri is positive definite (1151) implies that this equaling zero necessitates c j = 0 for all 
j = 1,..., N. Thus g only consists of the boundary terms, i.e. 

M 

g = ^2^div(-,yi)d u 

i=i 

from which one can show similarly that 

HsllAt # (R d ) ~ ^ *&div(yh ym)d m , 

l,m 

and since $div is also positive definite we must have d/ = 0 for all Z = 1,..., M. This completes the 
proof. □ 

Once (l23l) is solved, the resulting approximation decomposes as follows: 

N M N 

Sf — ^ Tr/??; (', dlj jCy -)- Qdivi.': yj )Hyj dj ^ ^ rljr j (•. Xj )Cj . 

j —1 j— 1 j— 1 

V V 

Pdiv^f Pcurl&\ 
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As a bonus, we get a stream function ip s t and velocity potential q s t satisfying 

4 = curl(i/> s t) + Vg s t. (24) 

Indeed, the identities m imply that such potentials are given by 

N M N 

^ S ‘ := -5Z curl W'> a: i) c j) -^2 curl (^(-’ x j) n yj) d j and H ■= -^2 V T ((j>(-,Xj)cj). 

3= 1 3=1 3=1 

4.2 Kernel Approximation with Curl-free Boundary Conditons 

We now focus on how to obtain a kernel-based approximation to the decomposition in Proposition 
[2j whose gradient term Vp is normal to the boundary. As in the previous section, we enforce full 
interpolation on a node set X and apply boundary conditions on a node set Y. The boundary 
conditions are imposed in this case by first projecting a kernel approximation s" onto the subspace 
of curl-free functions, and then setting all tangential components to zero pointwise. In d = 2 
dimensions, this is given by ty.P cur isf(yj) = 0 for all ijj £ Y, where t Vj is tangent to T at yj. As 
before, the Riesz representers give the basis functions one should consider: for full interpolation 
they are the same as the previous section, and the boundary-centered basis functions are of the form 
§curi{.-,yj)ty r Thus the interpolant is written as 


N M 

Sf = 'y ] 4>(-, Xj )c j + y ( &curl ('iVj )t yj dj . 
3=1 3 = 1 


(25) 


In the d = 3 case the two dimensional boundary leads to two basis functions at each shift on the 
boundary. For notational simplicity, we will continue with the d = 2 case here. 

The interpolation constraints give rise to a linear system similar to (1231) for determining the 


coefficients c j and dj: 


A B 
B t C 


c 

d 


where A is the matrix given in (1141) , B is given by 

*&curl (x± , yi )ty x 


B = 


$curl(x N ,yi)t 


Vl 


fix 

0 


^curl{.Xi,y M )t 


(26) 


UM 


*&curl (xN > yM )t 


VM 


and C is the M x M matrix with Cij = t J. & CU ri (jg, yj )t Vj .It can be shown using an argument similar 
to that in Lemma [2] that the linear functionals involved are linearly independent, which guarantees 
that the matrix in (1261) is symmetric and positive definite. The decomposition of the resulting kernel 
approximation is given by: 

N N M 

Sf — y ^ Xj ) Cj T y ( ^curz(*i Xj ) Cj T y ( &curl(’> yj)t Vj dj . 

3=1 3 = 1 3 = 1 


Pdiv Sf 


Pcurl S f 


In Section [5.21 we will show that Pdiv s" and Pcurisf approximate the terms from Proposition [5] 
Also one can use the form of the kernels (1161) to access potential functions i/> s n and q s n. 
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5 Error Estimates 


Our analysis follows the paradigm of RBF error estimates developed in recent years, where bounds 
on Sobolev functions having many zeros (the so-called “zeros lemmas,” or “sampling inequalities”) 
play a prominent role [23] . We will review the specific results we require below, and extend them 
slightly to suit our purposes. Next, we derive the error estimates in Sections [5.21 and [5731 

5.1 Zeros Lemmas 

The zeros lemmas involve bounding the norm of Sobolev functions that vanish on a set X = 
{xi ,..., Xn} C LI C in terms of the density of X in LI, which is quanitfied by the mesh norm: 

Hq := sup dist(x, X). 

The following is from [23] . with improvements in EH Theorem 4.6]. 

Proposition 6 . Let Ll C be a bounded domain with Lipschitz boundary. Let s £ R with s > d/2, 
and let p, £ K satisfy 0 < p < s. Also, let X C LI be a discrete set with mesh norm ho, sufficiently 
small. Then there is a constant depending only on LI such that if hn < Cq, and ifu£ H s (Ll) satisfies 
u\x = 0, then 

IMI/nqn) < Ch s n ( 27 ) 

where the constant C is independent of hn and u. 

This result can also be extended to manifolds in a straightforward way (see p3J Lemma 10]). Thus, 
if u £ H s ( T) satisfies u\y = 0, for 0 < p < s one has 

IM|.H>(r) — /i |l M llff s (r)- ( 28 ) 

Here the mesh norm hr for a finite set Y C T, is defined just as in the Euclidean case, the only 
difference being that distances are measured on the surface T. 

Note that the proposition above, the smoothness in the norm on the right-hand-side of the 
estimate is assumed to be high-enough so that the associated space of functions is continuous. 
However, such estimates hold for continuous functions in rougher norms, that is, if s > max{d/2,1} 
and u € H s (Ll) satisfies u\x = 0, ther0 

IMU 2 (fi) ^ Chn\u\ H iw 

If the underlying domain is a surface, by applying this estimate on patches, we get the following for 
continuous functions u : T —> R with zeros on Y CT: 

IMU 2 (r) < Chr\u\ H nry (29) 

Lastly, in what follows we will need zeros estimates in negative-indexed Sobolev norms. Note 
that if u € C(T) fl 7L 1 (T), then obviously u 6 H 1 (T) C 1L _1 (T). Thus we get 

IMU-i(r) = sup {u,ip) = |M|| 2(r) /IMIffi(r), 

IMI H i(r ) =1 

where since u € H 1 ( T) the supremum is achieved by choosing ip = u/||w||jji(r)- Thus if u vanishes 
on y, then with (l29l) we obtain 

IMU-qr) = ll u lli 2 (r)/ll u llffi(r) < C'ftr||«|U 2 (r)- (30) 

5 The proof of Proposition [§] involves local polynomial approximations on patches - in this case the polnomials are 
simply constants, which greatly simplifies the arguments. 
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5.2 Convergence with Divergence-free Boundary Conditions 

For the rest of the paper we assume that the RBF (f> is such that A/$(R d ) = H r (fl) with equivalent 
norms, the boundary T is smooth (at least C™’ 1 with 0 < r < m ), and that the mesh norms for the 
node sets X and Y (h^ and hr) are sufficiently small for the zeros lemmas to be applied. Further, 
we assume that g satisfies the condition (g, l) ri =0on each connected component of F. We begin 
with a basic interpolation estimate. 

Lemma 3. Let g satisfy 0 < g < r. Let s£ be the kernel approximation discussed in Section \4-l\ for 
a given f and g. Then for all f £ H r (f2) and g £ H T ~ 1 / 2 (r) we hav 

||f - Sf || H m(o) < Chf} (||f ||h t (Q) + llffll// T - 1 /2(r)) • 

Proof. Since f — s£ has zeros on X , we may apply Proposition [G] to get 

||f - SfllHM(a) < CHq M ||f - Sf|| H -(n). 

Now we use the extension operator. Since Ei\n = f and (PdivEf)\n = w, where w satisfies w n = g, 
then the data in the system used to determine s|, f (see ([23])) is the same as that of sj. Thus we 
get s£ = s|. f . This with (HHI) . the fact that H r (R d ) is norm equivalent to A/$(R d ), and and the 
continuity of E gives 

||f ~ 4llH’-(a) = \\Ei — s|; f ||h t (Q) < ll^f ~ s £f llH T (R< i ) — C\\Ef - s^ f || A q,(R<J) 

< < C||F/f< C (||f|| h^(£2 ) + ll^ll• 

This completes the proof. □ 

We continue our our analysis by showing that Pdi V s£ • n — g is small on the boundary. 

Lemma 4. Let g satisfy 0 < g < r. For all f £ H r (fl) and g £ H T ~ 1 ^ 2 (T) we have 

11 Pdiv Sf ■ n — g 11 nn 1 / 2 < Chf ^ (|| f ||(£2) + llfl’lli?’—V 2 ^)) • 

Proof. First assume that g > 1/2. Recall that Pdi V s\ ■ n = g on the node set Y C T by construction. 
Since the normals are assumed smooth and g — 1/2 > 0, we can apply (12811 to get 

llPd^sJ; • n - g\\ H „-i/ 2 ( r) < Chl~^~ 1/2 \\P d ivSf ■ n-gWHT-i/ 2 ^ 

< Ch T r A ' (|| Pdiv^ || H T_ !/ 2 (r) + llffllfl’’—V 2 ^)) • 


Applying the Trace Theorem and the fact that the H r (R“) norm bounds the H r (R“) norm gives us 
||-fdiu s f llH T - 1 / 2 (r) A C||P(jiii s f llH' r (n) ||Pdi!; s f IlfjT^gdj = C|| Pdiv s Ef 1 1 H T (R d ) — C ||&Ef IIH T (R d ) ’ 


where in the last two steps we used the fact that s£ = s^, f and that Pdi V is a projection on H T (R d ). 
The continuous embedding of A/*(R d ) into H T (R d ), the bounds (j2lj) . and continuity of E gives us 

ll s .Ef IIht(r<1) < C'||s| !f || J v* ( R- ) < C||Pf||^ (Rd ) < CIIMIIg T (R d ) — C||f|| H T (P)- 


This gives us the correct approximation orders down to g = 1/2. To get the estimates for 0 < 
g < 1/2, wc will measure the error in the H~ 1 ( T) norm, and then obtain the desired bound by 
interpolation. 

,: ‘Here and throughout, C is a constant independent of f, g, and the node sets. 
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Let V = {v £ H T 1 / 2 (r) : (v, l)|r 4 = 0 , 0 < i < K}, and note that this space is closed in the 
H t ~ 1 ' 2 (T) norm. Next consider the Banach space B := H T (fi) x V with obvious norm ||(f,g)||g := 
IIfIIH’-(n) + llffllif’— 1 / 2 (r)- Now define the linear map T : B —>• L 2 (T) given by T(f,g) := PdivS^-n-g. 
The argument above shows that 

||T||^ i2(r) < C7ip -1/2 . 

Similarly, considering T as a map from B to H~ l {Y ), the zeros estimate (1301) applies to the same 
arguments above to yield 

l|7 1 ||g_ s .,H-i(r) < C'/ip +1//2 . 

Estimates for the space iJ^ _ 1 / 2 (r) now follow from interpolation theory. Specifically, the identity 
for interpolation spaces in (3]) with 0 = 1/2 — /i gives us that 

Interpolation of operators (see, for example [J1 Proposition 14.1.5]) tells us that T maps B into 
H»- i/ 2 ( r ) with norm: 

imi B _>*,-i/ a( r) < < Ch r~^ 

This finishes the proof. □ 

Next, apply Proposition |T] to obtain sj; = w s t + Vp s t. Next we show that PdivSf approximates 

W„t. 

Lemma 5. Let 0 < /i < t. For all f £ H T (fl) and g £ H T ~ 1 ^ 2 (T) we have 

H-PdiuSf — w s t|| H „(Q) = ||P CMr iSf — Vp s t|| H n(Q) < Chi M (||f||H T (n) + ||fl r ||ir- r - 1 / 2 (r)) ■ 

Proof. The first equality follows easily from fact that Pdi V s\ — w s t = Vp s t — P cur is\. For the rest, 
note that P cu ws£ = Vq s t, where q s t is from (1241) . It follows that 

PdivSf = w s t + V(p.t — q s t), 

which is the decomposition in Proposition Q] applied to the function f = Pdi V s £. Letting v := 
Pdiv Sf — w s t, by dHJ) we get the bound 

ll v llL 2 (n) = ||V(p B t - g s t)|| L 2 (n) < C||Pdij; s f ■ n- ff|| ff -i/ 2 (r ). 

An application of Lemma [I] finishes the proof for the p = 0 case. For p, > 1, we can use dj) to get 

IMIsnqn) ~ lll v llln = II v IIl 2 (q) + ll v • n llIfM—!/2 (r) , 

where we used the fact that v is divergence-free and curl-free. After applying the bound on ||v|| L2 (q) 
above and the fact that v n = Pdi V s\ ■ n — g, we get 

ll v llH^(n) < C\\Pdiv s f ■ n — g||ffM-i/2(r)- 

Another application of Lemma [4] finishes the proof for 1 < p < r. The 0 < p < 1 case can be 
handled by interpolating the operator T between the ranges Lo(S^) and H 1 )!!), where T is given by 
T{ f, g) := Pdiv4 - w s t for (f, g) £ B. □ 

Now we are ready to prove one of our main results. 

Theorem 1. Let 0 < p < r. Given f £ H T (fl) and admissible g £ H T ~ 1 ' 2 (T), we denote the 
decomposition off from Proposition Q] as f = Wf + Vpf. Then we have 

ll-PdwSf — w f [| H /*(n) < C (hi M + hp A ') (||f||H T (n) + IlsHip— i/ 2 (r)) • 
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Proof. We begin with a triangle inequality and an application of Lemma [5] 

H-PdruSf - Wf ||HM(fi) < ||w s t - Wf|| H ^(0) + C/lp M ( ||f ||h-(O) + ||s||H’- 1 /2(r)) • 

Next we bound ||w s t — Wf|| H /qn)- Note that s£ — f = (w s t — Wf) + V(p s t — p?) decomposes s£ — f 
as in Proposition Q] with g = 0. Applying Proposition [4] to f — s£, we get that w s t — Wf = curl(-0) 
with xp satisfying m, which yields: 

l|w s t - w f || HM(n) = ||curl(i/>)|| H (‘+i(n) < C||4 - f|| H M(n)- 

An application of Lemma [3] finishes the proof. □ 

Since P CU ris\ — Vpf = s£ — f + Wf — Pdi V s\ , similar estimates hold for the curl-free part. 


5.3 Convergence with Curl-free Boundary Conditions 

Now we focus on the decomposition in Proposition [2] Recall that there is a projector P n that 
projects f onto the curl-free term in this decomposition, and that sj? denotes the kernel interpolant 
from Section l4~2l whose tangential components of P cllr ;sf? are forced to vanish on the node set hcP. 
Showing that P cur /sj? approximates P n f uses arguments similar to those in the preceeding section, 
thus we provide only the aspects of the proof that are significantly different 

First, we have a lemma, whose proof we omit since the arguments are similar to those of Lemma 
[3]- the most major difference here is that the proof requires an extension E so that sf? = s£. f , and 
such an extension exists by Lemma 1 and the remark proceeding it. 

Lemma 6. Let /i satisfy 0 < p < r. Then for all f £ H T (fl) we have 

||f - Sf |!h^(Q) < Ch T n ^ ||f|| H ' r (n)• 

Next we have a lemma analogous to Lemma [5j 

Lemma 7. Let 0 < p, < t. Then for all f £ H T (fl) we have 


|P n s" - P CU riSf ||hm(Q) < Ch T r 


iH’-(fi)- 


Proof. We will use the tangential trace operator y t , which is defined on smooth vector fields as 
7 t v := v|r x n. By [TS] Theorem 2.11, page 34], this extends to a continuous map defined on L 2 (f2) 
vector fields with bounded curl (in L 2 ) to the space H -1 / 2 (r), and the following Green’s formula 
holds: 

(curl(v),g) - (v,curl(g)) = ( 7 t v,g) Vg £ H 1 (fl). (31) 

The first step is to transfer the problem to the boundary by showing that 

||Pn s f — Pcurl Sf 11H/ 2 (H) ^ CHytPctiriSf II h ^- 1 / 2 (r) ■ (32) 


For brevity, we let v = P n sf? — P cu wsf?. The identity 

Pn Sf - PcurlS f = Pdiv Sf - P^ Sf 

implies that v £ curl(H 1 (fl)), so by Proposition [3] v has a potential xp satisfying 


II'0Hhi(£ 2) < C||v|| L 2 (n). 

With this, we can apply (1311) with g = xp to get the inequality 

IMlLlfi) = l(7tv,-0)| < ||7tv|| H -i/2(r)ll^|lHi/2(r) 

< ll 7 tv|| H -i/ 2 (r)ll' i / , llH 1 (n) < C'||7tv|| H -i/2( r )||v|| L2 ( r ). 
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Since 7 tP n s" = 0 , we obtain (13^1) when p = 0: 

IMlL 2 (fi) < C\\j t PcuriSf || H -i/2( r ). (33) 

For p, > 1, we use ( 0 . Using (l33l) and the fact that v is both divergence-free and curl-free, we get 

11 v II (O) — f~lll u lllt = C(ll v llL 2 (n) + II Ttvll hm —!/2 (r)) — ^ll7t v llHf'- 1 / 2 (r) = C|| 7 t-Pcur;s f || Ht ,_i/ 2 ^ r j. 

This proves (l32l) for p = 0 and 1 < p < r. By design 7 t P cur isf has many zeros on T, which makes 
this situation very similar to that in Lemma [H whose arguments can be repeated to arrive at the 
bound 

||7t-fom-zSf ||h^—!/ 2 (r) < C^-r M II^II(O)- 

The case 0 < /i < 1 can now be handled by operator interpolation. This finishes the proof. □ 

With these results, one can now construct an argument very similar to the proof of Theorem [T] 
to arrive at the result below, which we state without proof. 

Theorem 2. Let 0 < p < t. Then for all f £ H T (fl) we have 

ll^nf - P™ws f n || H ,(fi) < c (h^ + h T r ~») ||f|| H r(n). 

Remark 1. We heavily relied on the fact that given f £ H r (f2), we are guaranteed potential functions 
having the appropriate smoothness (assuming T is smooth enough). We are not aware of such a result 
for functions in native spaces associated with C°° kernels, even for very smooth domains. However, 
convergence results for the decompositions treated here can be derived for C°° kernels, assuming that 
all potentials (or their components) reside within J\fwhere 4> = —A <j>. 


6 Numerical Examples 

In this section we illustrate the methods described previously with numerical experiments. We start 
with the following target function: 

f = curl(cos(2(x 2 + y 2 ))) + Vp, (34) 


where p is the MATLAB peaks function, and consider f on the annulus f 1 centered at the origin with 
inner radius .75 and outer radius 2 (see Figure 1(a) I. This function on U has the property that the 


Leray projection, P^f, is equal to curl(cos(2(ar+p 2 ))), and in what follows we will compare P/ f to 
PdivS f ■ We used the freely available distmesh package to generate quasi-uniformly spaced nodes on 
Q [24] for the experiments. Eight nodes sets were generated with the number of full-interpolation 
centers ranging from N = 615 to N = 11210, and the number boundary centers ranging in cardinality 
from M = 115 to M = 521. An example node set with N = 1276 is pictured in Figure [l(b) | In 
every experiment, we enforced full-interpolation at all centers, including the boundary sites. MATLAB 
files containing the nodes used and other useful files can be downloaded from m ■ To generate our 
matrix-valued kernels, we used the scalar Matern kernel </> given by 


<t>(r) = -—e~ r (r 5 + 15 r 4 + 105r 3 + 420r 2 + 945r + 945), 
945 


where r = r(x,y) = e \/x 2 + y 2 . The free parameter e, known as the shape parameter, affects the 
stability and accuracy of the method. The shape parameter remained fixed at e = 5 throughout 
our experiments, which kept the computations relatively stable. The two dimensional version of 
this kernel, <f>(\Jx 1 + y 2 ), satisfies </>(w) = C(1 + |cc| 2 ) -13 / 2 , where C is a constant, which means in 
particular that the matrix kernel 4> satisfies (BUI) with r = 5.5. 
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Figure 1: The domain and target field f used in the first experiment. 


We measured the relative error ||P^Sf — P^Wi^x)/\\P l %where X is the finest node set 
of those described above (i.e. with j/X = 11210) and the norm is given by 


\S\k 2 {x) = 


\ 


1 

w 


E 

Xi£X 


g(2b)| 2 - 


The error between the generalized interpolant s£ and f was recorded similarly. Lemma [3] and 
Theorem [T] dictate that the L, 2 (P) errors should all decay like (D(h 5 5 ). Since our nodes are very 
uniform, || • ||^ 2 (.y) ~ II ■ ||L 2 (n)> so observing 0(h 5 ' 5 ) would confirm these results. Due to the quasi¬ 
uniformity of the nodes, the mesh norm ft behaves asymptotically like 1/y/N, where N is the number 
of nodes in a given node set. A loglog plot of error versus 1/y/N is given in Figure 3(a)| where it 
can be seen that the error for the Leray projection appears to converge slightly faster than 0(h 5 - 5 ). 


In the next experiment, we computed the full Helmholtz-Hodge decomposition (HHD) of f on a 
slightly more complicated domain, and in the process obtained evidence for the bound in Theorem 
[2j Recall that the full HHD is given by 


f = P„f + P L f + Vft, 


(35) 


where P n f is the curl-free normal component of f from Proposition (2j Pz,f is the Leray projection, 
and ft, is a harmonic function. We used the same target function (1341) . but on the domain pictured 
in Figure 4(a) As in the previous test, several quasi-uniform node sets were generated using the 


distmesh package with sizes ranging from N = 486 to N = 16882 (see [16]). Samples of f at these 
sites were used to obtain approximations to each term in (l35l) using the method described below. 

The first step of the two-step process is to construct an interpolant of f with curl-free boundary 
conditions of form (ESI) that solves the system ESI). Let s" denote this interpolant and note that 
Pcurisf approximates P n f- Second, decompose Pdiv s" to approximate Pl f and Vft by using an 
interpolant with divergence-free boundary conditions of the form (1221) that solves (E3l) (with g = 0 
and f replaced by Pdiv sj?). Denote this interpolant by s£, and note that Pdiv Sf ~ P^f and P CU riSf ~ 
Vft. These steps give approximations to the three components of the decomposition of f, which are 
plotted in Figure 01 together with contour plots of the corresponding potential functions. 

With regard to convergence, we did not measure the error directly because the exact decompo¬ 
sition for f on this domain is unknown to us. Nevertheless, we estimated the rate of convergence 
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Figure 2: The kernel decomposition of f using s}; = Pdi V Sf + PcurlSf- The contours represent the 
potentials ip s t and q s t. 
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Figure 3: Convergence results for each numerical experiment. The vertical axis gives the logarithm 
of the relative ^(X) error (base 10), and the horizontal axis gives TV on a log 10 scale. 



Figure 4: The kernel approximation of the full HHD for the target field f (l34l) . with contours of each 
term’s scalar potential. 
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by using each approximation on the finest node set as proxies for the true solution. To measure the 
error corresponding to P n f, for example, we used ||P CM riSf — Vp||f 2 (x) where Vp is the kernel ap¬ 
proximation to P n f on the finest node set X (with ffX = 16882). We also tested the error between 
the generalized interpolant sj? and f. Lemma [3] and Theorem [2] dictate that the L 2 (fl) errors should 
all decay like 0(h 55 ). A loglog plot of error versus 1/y/N ~ h is given in Figure [3(b)| where the 
errors seem to be converging like 0(h 5 ' 5 ). 

7 Concluding Remarks 

Decompositions with other boundary conditions are certainly also possible. If no boundary con¬ 
ditions are specified, one can find an interpolant Sf using only shifts of positive definite kernel 
$ = —A cj)I. Enforcing Sf|x = f|x leads to a positive definite system, and since $ = + & C urh 

Sf decomposes trivially. This idea was used in a decomposition technique using thin plate splines 
introduced in earlier work [lj. For other boundary conditions, if the functionals associated with 
the interpolation and boundary conditions are linearly independent and the Reisz representers are 
chosen as basis functions, then the kernel decomposition can be constructed. In this way, one could 
impose a whole host of boundary conditions in vector decomposition problems, and do so in a natural 
way. 


References 

[1] L. Amodei and M. N. Benbourhim. A vector spline approximation. J. Approx. Theory , 67(1):51- 
79, 1991. 

[2] R. K. Beatson, J. Levesley, and C. T. Mouat. Better bases for radial basis function interpolation 
problems. J. Comput. Appl. Math., 236(4):434-446, 2011. 

[3] M. Benzi, G. H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta 
Numerica, 14:1-137, 2005. 

[4] S. C. Brenner and L. R. Scott. The mathematical theory of finite element methods, volume 15 
of Texts in Applied Mathematics. Springer-Verlag, New York, second edition, 2002. 

[5] M. D. Buhmann. Radial basis functions: theory and implementations, volume 12 of Cam¬ 
bridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 
Cambridge, 2003. 

[6] A. J. Chorin. Numerical solutions of the Navier-Stokes equations. Math. Comput., 22:745-762, 
1968. 

[7] R. Dautray and J.-L. Lions. Mathematical analysis and numerical methods for science and 
technology. Vol. 3. Springer-Verlag, Berlin, 1990. Spectral theory and applications, With the 
collaboration of Michel Artola and Michel Cessenat, Translated from the French by John C. 
Amson. 

[8] F. M. Denaro. On the application of the Helmholtz-Hodge decomposition in projection methods 
for incompressible flows with general boundary conditions. International Journal for Numerical 
Methods in Fluids, 43(l):43-69, 2003. 

[9] E. Deriaz and V. Perrier. Orthogonal Helmholtz decomposition in arbitrary dimension using 
divergence-free and curl-free wavelets. Appl. Comput. Harmon. Anal., 26(2):249-269, 2009. 


19 




[10] P. Farrell and H. Wendland. RBF multiscale collocation for second order elliptic boundary 
value problems. SIAM J. Numer. Anal., 51(4):2403-2425, 2013. 

[11] G. E. Fasshauer. Meshfree approximation methods with MATLAB, volume 6 of Interdisciplinary 
Mathematical Sciences. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2007. With 
1 CD-ROM (Windows, Macintosh and UNIX). 

[12] W. Freeden and T. Gervens. Vector spherical spline interpolation—basic theory and computa¬ 
tional aspects. Math. Methods Appl. Sci., 16(3):151 183, 1993. 

[13] E. Fuselier, T. Hangelbroek, F. J. Narcowich, J. D. Ward, and G. B. Wright. Localized bases 
for kernel spaces on the unit sphere. SIAM J. Numer. Anal, 51 (5):2538—2562, 2013. 

[14] E. Fuselier and G. B. Wright. Scattered data interpolation on embedded submanifolds with 
restricted positive definite kernels: Sobolev error estimates. SIAM J. Numer. Anal., 50(3): 1753— 
1776, 2012. 

[15] E. J. Fuselier. Improved stability estimates and a characterization of the native space for 
matrix-valued RBFs. Adv. Comput. Math., 29(3):269-290, 2008. 

[16] E. J. Fuselier. Matlab examples for ”A radial basis function method for computing helmholtz- 
hodge decompositions”. http://math.highpoint.edu/~efuselier/RBFVectorDecomposition/, 
page (accessed 02/24/2015), 2015. 

[17] E. J. Fuselier and G. B. Wright. Stability and error estimates for vector field interpolation and 
decomposition on the sphere with RBFs. SIAM J. Numer. Anal., 47(5):3213-3239, 2009. 

[18] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations, volume 5 
of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. Theory and 
algorithms. 

[19] P. Grisvard. Elliptic problems in nonsmooth domains, volume 24 of Monographs and Studies in 
Mathematics. Pitman (Advanced Publishing Program), Boston, MA, 1985. 

[20] D. Handscomb. Local recovery of a solenoidal vector field by an extension of the thin-plate 
spline technique. Numer. Algorithms , 5(1-4):121-129, 1993. Algorithms for approximation, III 
(Oxford, 1992). 

[21] J.-L. Lions and E. Magenes. Non-homogeneous boundary value problems and applications. Vol. I. 
Springer-Verlag, New York, 1972. Translated from the French by P. Kenneth, Die Grundlehren 
der mathematischen Wissenschaften, Band 181. 

[22] F. J. Narcowich and J. D. Ward. Generalized Hermite interpolation via matrix-valued condi¬ 
tionally positive definite functions. Math. Comp., 63(208):661-687, 1994. 

[23] F. J. Narcowich, J. D. Ward, and H. Wendland. Sobolev bounds on functions with scattered 
zeros, with applications to radial basis function surface fitting. Math. Comp., 74(250):743-763 
(electronic), 2005. 

[24] P.-O. Persson and G. Strang. A simple mesh generator in Matlab. SIAM Rev., 46(2):329-345 
(electronic), 2004. 

[25] K. Polthier and E. Preuss. Identifying vector fields singularities using a discrete Hodge decom¬ 
position. In H. C. Hege and K. Polthier, editors, Visualization and Mathematics III. Springer 
Verlag, Berlin, 2002. 


20 


[26] D. Schrader and H. Wendland. A high-order, analytically divergence-free discretization method 
for Darcy’s problem. Math. Comp., 80(273):263-277, 2011. 

[27] G. Schwarz. Hodge decomposition—a method for solving boundary value problems, volume 1607 
of Lecture Notes in Mathematics. Springer-Verlag, Berlin, 1995. 

[28] E. M. Stein. Singular integrals and differentiability properties of functions. Princeton Mathe¬ 
matical Series, No. 30. Princeton University Press, Princeton, N.J., 1970. 

[29] R. Temam. Sur l’approximation de la solution des equations de Navier-Stokes par la methode 
des pas fractionnaries, II. Arch. Ration. Mech. Anal., 33:377-385, 1969. 

[30] H. Wendland. Scattered data approximation, volume 17 of Cambridge Monographs on Applied 
and Computational Mathematics. Cambridge University Press, Cambridge, 2005. 

[31] H. Wendland. Divergence-free kernel methods for approximating the Stokes problem. SIAM J. 
Numer. Anal., 47(4):3158-3179, 2009. 


21 


